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ABSTRACT 



Context. The solar torsional oscillations, i.e., the perturbations of the angular velocity of rotation associated with the eleven-year 
activity cycle, are a manifestation of the interaction among the interior magnetic fields, amplified and modulated by the solar dynamo, 
and rotation, meridional flow and turbulent thermal transport. Therefore, they can be used, at least in principle, to put constraints on 
that interaction. Similar phenomena are expected to be observed in solar-like stars and can be modelled to shed light on analogous 
interactions in different environments. 

Aims. The source of the torsional oscillations is investigated by means of a model for the angular momentum transport within the 
convection zone. 

Methods. A description of the torsional oscillations is introduced, based on an analytical solution of the angular momentum equation 
in the mean-field approach. It provides information on the intensity and location of the torques producing the redistribution of the 
angular momentum within the convection zone of the Sun along the activity cycle. The method can be extended to solar-like stars for 
which some information on the time-dependence of the differential rotation is becoming available. 

Results. Illustrative applications to the Sun and solar-like stars are presented. Under the hypothesis that the solar torsional oscillations 
are due to the mean-field Lorentz force, an amplitude of the Maxwell stresses \B t B$\ ^ 8x 10 3 G 2 at a depth of ~ 0.85 R e at low latitude 
is estimated. Moreover, the phase relationship between B T and B^ can be estimated, suggesting that B r B,p > below ~ 0.85 R a and 
B r B(6 < above. 

Conclusions. Such preliminary results show the capability of the proposed approach to constrain the amplitude, phase and location 
of the perturbations leading to the observed torsional oscillations. 

Key words. Sun: rotation - Sun: activity - Sun: magnetic fields - Sun: interior - stars: rotation - stars: activity 
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1. Introduction 

Doppler measurements of the surface rotation of the Sun show 
bands of faster and slower zonal flows that appear at midlat- 
itudes and migrate toward the equator with the period of the 
eleven-year cycle, accompanying the bands of sunspot activity. 
The amplitude of such velocity perturbations, called torsional 
oscillations, is of ~ 5 m s _1 and faster rotation is observed on 
the side equatorward of the sunspot belt (Howard & LaBonte 
1980). Helioseismology has revealed that the torsional oscilla- 
tions are not at all a superficial phenomenon but involve much 
of the convection zone, as shown, for example, by Howe et al. 
(2000), Vorontsov et al. (2002), Basu & Antia (2003) and more 
recently by Howe et al. (2005, 2006). The amplitude of the an- 
gular velocity variation is 6Q./2n ~ 0.5 - 1 nHz at least down 
to 10%- 15% of the solar radius, although the precise depth of 
penetration of the oscillations is difficult to establish given the 
present uncertainties of the inversion methods in the lower half 
of the solar convection zone (e.g., Howe et al. 2006). In addition 
to such a low-latitude branch of the torsional oscillations, he- 
lioseismic studies have detected the presence of a high-latitude 
branch (above ~ 60° latitude) that propagates poleward and the 
amplitude of which is about 6Q./2tt ~ 1 — 2 nHz (see also, e.g., 
Toomre et al. 2000; Basu & Antia 2001). Such a branch seems to 



propagate almost all the way down to the base of the convection 
zone. 

A general description of the perturbation of the angular ve- 
locity of the torsional oscillations, given the present accuracy 
of the observations, is provided by the simple formula (cf., e.g., 
Vorontsov et al. 2002; Howe et al. 2005): 

cj(r, 6) = A (c) (r, 8) cos(o-f) + A (s) (r, 0) sin(crf) = 

= A(r,6)sm[o-t + <f>(r,0)], (1) 

where r is the distance from the centre of the Sun, 9 the co- 
latitude measured from the North pole, cr the frequency of the 
eleven-year cycle, t the time, and the amplitude functions A^ = 
A sin (p, A (s) = A cos <p depend on the amplitude A and the initial 
phase (f>. Moreover, the velocity perturbation is symmetric with 
respect to the equator: 



(o(r,0) =u>(r,n-ff). 



(2) 
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Several models have been proposed to interpret the torsional 
oscillations beginning with the pioneering work by Schussler 
(1981) and Yoshimura (1981) who considered the Lorentz force 
associated with the magnetic fields in the activity belts as the 
cause of the velocity perturbations observed in the solar photo- 
sphere. Later models, based on the effects of the Lorentz force on 
the turbulent Reynolds stresses, were proposed, by, e.g., Kiiker 
et al. (1996) and Kichatinov et al. (1999), following an original 
suggestion by Rudiger & Kichatinov (1990). Spruit (2003) pro- 
posed that the low-latitude branch of the torsional oscillations is 
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a geostrophic flow driven by temperature variations due to the 
enhanced radiative losses in the active region belts. 

More recent works by Covas et al. (2004, 2005) present mod- 
els based on the simultaneous solution of non-linear mean-field 
dynamo equations and the azimuthal component of the Navier- 
Stokes equation with a uniform turbulent viscosity. They repro- 
duce the gross features of the torsional oscillations and of the 
solar activity cycle with an appropriate tuning of the free param- 
eters. Rempel (2006, 2007) considers the role of the meridional 
component of the Navier-Stokes equation in mean-field models 
and finds that the perturbation of the meridional flow cannot be 
neglected in the interpretation of the torsional oscillations. His 
models suggest that the low-latitude branch of the torsional os- 
cillations cannot be explained solely by the effect of the mean- 
field Lorentz force, but that thermal perturbations in the active 
region belt and in the bulk of the convection zone do play an 
active role, as proposed by Spruit (2003). 

In the present study, the angular momentum conservation is 
considered and the relevant equation in the mean-field approxi- 
mation is solved analytically for the case of a turbulent viscosity 
that depends on the radial co-ordinate. A general solution is de- 
rived independently of any specific dynamo model, allowing us 
to put constraints on the localization of the torques producing the 
torsional oscillations. An illustrative application of the proposed 
methods is presented using the available data. 

The observations of young solar-like stars by means of tomo- 
graphic techniques based on high-resolution spectroscopy have 
recently provided evidence for time variation of their surface dif- 
ferential rotation (see, e. g., Donati et al. 2003; Jeffers et al. 
2007). Lanza (2006a) has recently shown how such variations 
can be related to the intensity of the magnetic torque produced 
by a non-linear dynamo in their convective envelopes, in the case 
of rapidly rotating stars for which the Taylor-Proudman theorem 
applies. In the near future, the possibility of measuring the time 
variation of the rotational splittings of p-mode oscillations in 
solar-like stars may provide us with information on the changes 
of their internal rotation, although with limited spatial resolu- 
tion. In the present study, we extend the considerations of Lanza 
(2006a) to the case of a generic internal rotation profile, not nec- 
essarily verifying the Taylor-Proudman theorem, to obtain hints 
on the amplitude of the torque leading to the rotation change. 

2. The model 

2.1. Hypotheses and basic equations 

We consider an inertial reference frame with the origin in the 
barycentre of the Sun and the z-axis in the direction of the 
rotation axis. A spherical polar co-ordinate system (r, 6, <p) is 
adopted, where r is the distance from the origin, 8 the co-latitude 
measured from the North pole and (p the azimuthal angle. We as- 
sume that all the variables are independent of if and that the solar 
density stratification is spherically symmetric. 

The equation for the angular momentum conservation in 
the mean-field approach reads (e.g., Riidiger 1989; Riidiger & 
Hollerbach 2004): 



d 

— (pr 2 sin 2 0Q) + V -0 = 0, 
ot 



(3) 



where p(r) is the density, £l(r, 6, t) the angular velocity and the 
angular momentum flux vector given by: 

= (pr 2 sin 2 6Q.)u (m) + 



y sin 

+r sin eipu'u'J - ——{BB V + (B'B'J), 



where M (m) is the meridional circulation, u' the fluctuating ve- 
locity field, p the magnetic permeability, B the mean magnetic 
field and B' the fluctuating magnetic field; angular brackets in- 
dicate the Reynolds average defining the mean-field quantities. 
The Reynolds stresses can be written as: 



, , / dui du;\ 

(p M(M .) = -, t + _ + A„ 



(5) 



where u is the mean flow field, rj t (r) is the turbulent viscosity, 
assumed to be a scalar function of r only, and A/y indicates the 
non-diffusive part of the Reynolds stresses due to the velocity 
correlations in a rotating star (see Riidiger 1989; Riidiger & 
Hollerbach 2004, for details). The conservation of the total an- 
gular momentum of the convection zone implies: 



® r = for r = fb, Rq, 



(6) 



where r\, is the radius at the lower boundary of the convection 
zone and R Q is the radius of the Sun. 

The equation for the conservation of the angular momentum 
can be recast in the form: 



<9Q 
~dt 



1 d 



, dQ. 

■ n, 

dr 



pr 2 (1 - p 2 ) dp 



2\2 



(i -aO 



3D. 
dp 



S, 



where p = cos 8 and the source term S is given by: 
V-r 



pr 2 (l - p 2 Y 
and r is a vector whose components are: 



(7) 



(8) 



T; = rsinf 



a,v - I (biB v + (b;b;>) 

A' 



+ pr 2 sin 2 6>nw (m )i. (9) 



The boundary conditions given by Eq. (6) can be written as: 

^ = for r = r b , R e , (10) 
or 

when we assume r r = at the surface. Note that helioseismic 
measurements indicate the presence of a subsurface shear layer 
with f: < at low latitudes (Corbard & Thompson 2005), but 
we prefer to adopt the stress-free boundary condition (10) at the 
surface to ensure the conservation of the total angular momen- 
tum of the convection zone in our model. 

The solar angular velocity can be split into a time- 
independent component Qo and a time-dependent component cj, 
i.e., the torsional oscillations: 



£l(r,p, f) = £lo(r,p) + io(r,p, t). 

The equation for the torsional oscillations thus becomes: 
du)\ 



(11) 



dco 1 d 
m 1 



pr 2 (1 - p 2 ) dp 



op 



= S U 



(4) Si = 



where the perturbation of the source term is: 

V-Ti 



pr 2 (\-p 2 Y 



(12) 



(13) 
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with 



T\i = rsin# 



where A* and m^ ); are the time-dependent perturbations of the 
non-diffusive Reynolds stresses and of the meridional circula- 
tion, respectively, and fio is the average of the solar angular ve- 
locity over the convection zone. Note that the Maxwell stresses 
appear in Eq. (14), but not in the corresponding equation for Qo 
because the solar magnetic field has no time-independent com- 
ponent. Moreover, in deriving Eq. (14), the variation of the an- 
gular velocity over the convection zone has been neglected in 
the term containing the perturbation of the meridional circula- 
tion since M <sc Q (cf., e.g., Rudiger 1989; Rempel 2007). 
Eq. (12) must be solved together with the boundary conditions: 



A (p) 



i(B^ + CB;s;>) 



+ pr 2 sin 2 0Q o <V (14) 



,(p) 



eigenvalues verify the inequality: A n o < A n \ < ...A„k < A nk+ \ < ... 
and that the eigenfunction £ nk has k nodes in the interval [r b ,R B ] 
for each n. For n = 0, the first eigenvalue corresponding to the 
eigenfunction foo is zero an d the eigenfunction vanishes at all 
points in [r b ,R & ], as it is evident by integrating both sides of 
Eq. (19) in the same interval, applying the boundary conditions 
(20) and considering that £oo has no nodes. For n > 0, all the 
eigenvalues A„q are positive, as can be derived by integrating 
both sides of Eq. (19) in the interval [rb,R Q ], applying the 
boundary conditions (20) and considering that £„o has no nodes. 
In view of the inequality given above, all the eigenvalues A nk are 
then positive for n > 0. Moreover, it is possible to prove that 
A„'k > A n t \fn'>n and that (see Smirnov 1964b): 



{£f)p + n(n + 3)q (*f) P + n(n + 3)6 
— < A nk < 



M 



m 



(21) 



da) 
~dr~ 



= for r = r h , R . 



(15) 



2.2. Solution of the angular momentum equation 



The general solution of Eq. (12) with the boundary conditions 
(15) can be obtained by the method of separation of the variables 
and expressed as a series of the form (cf., e.g., Lanza 2006b, and 
references therein): 



where P, Q and M are the maximum values of the functions 
r A r] t , r 2 77 t and pr 4 in the interval [r b ,Ro\, respectively, and p, q 
and m their minimum values in the same interval, respectively; 

and / = f ° ^Jp/r/tdr. For k » 1 the asymptotic expression for 
the eigenfunction £ n k is (cf., e.g., Morse & Feshbach 1953): 




(22) 



co(r,p,t) = ]T Y^ankitXnk^P^Kp), 



(16) 



The time dependence of the solution (16) is specified by the 
functions a„k that are given by: 



«=0,2,4,-- k=0 



da nk 



where a„k(t) and £„k(r) are functions that will be specified below ^ f + ^nk&nk(f) - Pnkif), 



(23) 



and P^' ) ( J u) are Jacobian polynomials, i.e., the finite solutions 
of the equation: 



where the functions fi„k appear in the development of the pertur- 
bation term S \ : 



d_ 



2x2' 



dp 



+ n(n + 3)(l-// 2 )P< U) = 0, (17) S r(r,pj) = J]^ nk ^ nk (r)^(p.), 



(24) 



in the interval -1 < p < 1 including its ends (cf., e.g., Smirnov 
1964a). The Jacobian polynomials form a complete and orthog- 
onal set in the interval [—1,1] with respect to the weight function 
(1 -p 2 ). Only the polynomials of even degree appear in Eq. (16) 
because the angular velocity perturbation is symmetric with re- 
spect to the equator (see Eq. 2). For n » 1 the asymptotic ex- 
pression of the Jacobian polynomials is (see, e.g., Gradshteyn 
& Ryzhik 1994): 



« k 

and are given by (cf. Lanza 2006b): 
_ (In + 3)(» + 2) 

Pnk — ~ — X 



8(n + 1) 



P<, u) (cos0) = 



cos [(« + \)6 - f ] 
V7rn[sin(f)cos(f)] 3/2 



X X Xi P, " 4(1 ~ ^ 2)S l(r '^ K^^Wdrdp, (25) 
The solution of Eq. (23) is: 



+ 0(n~) 



(18) 



a nk (i) = a nk (0) + exp(-i, 



nkt) f Pnk(t' 

Jo 



)exp(A nk t')dt'. (26) 



The functions ( nk are the solutions of the Sturm-Liouville 
problem defined in the interval r/, < r < R e by the equation: 



d 
dr 



r T] V 



d£ 



nk 

d7 



- n(n + 3)r rj^nk + Kkpr Cnk = 



with the boundary conditions (following from Eq. 15): 

— — = Oatr = r b ,R Q . 
dr 



(19) 



(20) 



We shall consider normalized eigenfunctions, i.e.: 
° pr 4 £ 2 k dr = 1. The eigenfunctions £ nk for a fixed n, 
form a complete and orthonormal set in the interval [r\,, R Q ] with 
respect to the weight function pr 4 that does not depend on n. We 
recall from the theory of the Sturm-Liouville problem that the 



which allows us to specify the general solution of Eq. (12) 
with the boundary conditions (15) when the perturbation term 
S i(r,p, t) and the initial conditions are given. 

2.3. Solution for the solar torsional oscillations 

To find the solution appropriate to the solar torsional oscillations 
as specified by Eqs. (1) and (2), it is useful to derive an alter- 
native expression for the functions f3 nk as follows. Substituting 
Eq. (13) into Eq. (25) and taking into account that the element 
of volume is dV = -2nr 2 drdp, the r.h.s. of Eq. (25) can be recast 
in the form of a volume integral extended to the solar convection 
zone: 



Pnk = F„ f(V • r^UP^dV, 
Jv 



(27) 
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where: 



F„ = 



(In + 3)(n + 2) 
I6n(n + 1) 



(28) 



is a factor coming from the normalization of the Jacobian poly- 
nomials. It is possible to simplify further Eq. (27) by considering 
the identity: 



(29) 



Integrating both sides of (29) over the volume of the convec- 
tion zone and considering that the integral of the l.h.s. vanishes 
thanks to the Gauss's theorem and the condition that the radial 
component of the stresses T\ r is zero on the boundaries, we find: 



(30) 



The time dependence of the solar torsional oscillations spec- 
ified in Eq. (1) suggests to consider a similar dependence for the 
perturbation term: 



T\(r,n, f) = T, (r,//) cos(crf) + t, (r,fi) sin(crf), 



(31) 



from which a similar expression for the follows by substi- 
tution into Eq. (27). If we put such an expression for /3 n k into 
Eq. (26) and perform the integrations with respect to the time, 
we find the stationary solution: 



a„k(t) 



F„ 



^UP ( n' l) {^-rf+crV-Tf)dV 



COs(crf)+ 
sin(crf) 



(32) 



This expression can be substituted into Eq. (16) to give the angu- 
lar velocity perturbation. It can be written in the form of Eq. (1) 
with: 

A (e) 0,//) = j [Gi(r,//, /,//)V • rf - G 2 (r,fx, r',//)V • t\ s) ] dV, 
A {s \r,n) = J [G x (r,n, r',/i')V ■ t\ s) + G 2 (r,n, r',fi')V ■ rf] dV, 



Eqs. (33) can be used to compute the angular velocity per- 
turbation when n is known. Note that in the case in which the 
Lorentz force due to the mean field and the meridional flow are 
the only sources of angular momentum redistribution, Eq. (14) 
gives: 



V • n = -^Bp • V(r sin0B ) + 2pClo(r sin0) M(m)s , 



(35) 



where B p is the mean poloidal magnetic field and M( m ) S the com- 
ponent of the meridional flow in the direction orthogonal to the 
rotation axis. To obtain Eq. (35) we made use of the solenoidal 
nature of the mean poloidal field and of the continuity equation 
for the meridional flow. 

2.4. Localization of the source of the torsional oscillations in 
the solar convection zone 

The results derived above allow us to introduce methods to lo- 
calize the torques producing the torsional oscillations in the con- 
vection zone. Suppose that the observations provide us with the 
functions A (c) (r,yu) and A (s) (r,//) appearing in Eq. (1). The func- 
tions a„k(t) can be written as: 



a^l cos(crf) + a { °l sin(crf), 



(36) 



where the constants a^ s) are given by: 



a 



(c,s) _ 
nk 



= 2nF n £ ° £ A M (r,n)pr\\ - /^P^W^u. (37) 



Similarly, we can write: 

Pnk(t) = b (C 2 COS(crf) + b ( *l Sin(crf), 

with the relationships: 



nk ' 

nk 



Anka,,,. + era 



A n kCl 



nk 
nk 



era 



nk ' 



nk ' 



(38) 



(39) 



that follow from Eq. (23). 



The divergence of t, c ' s) can be obtained from Eqs. (13) and 



(24) as: 



V • T (C ' S) 

(33) V T ' 



(40) 



where the symbol dV means that the volume integration is to 
be performed with respect to the variables r' and //, and the 
functions G\ and G2 are Green functions defined as: 



G l (r, i i,r',n') = Y j F n Y j 



&nk 



Moreover, it is possible to construct a localized estimate of the 
perturbation term T\ by considering a function f(r,fi) the gradi- 
ent of which is different from zero only within a given volume 
Vf. It can be developed in series of the eigenfunctions in the 
form: 



xU{r)P ( n' l) mnk{r')P ( n U) <4i') 

CO CO 

n k nk 

xU(r)P n l ' l Xki)U(r')P n ll \v')- 



f(r,n) = Y J YjC nk U(r)P n U) (M), 

n k 



(34) C nk 



The Green functions are continuous with respect to the argu- 
ments r, /j., r', /i', but their partial derivatives with respect to r', 
// have discontinuities of the first kind in the points where r — r' 
or [i = yu' '. The convergence of the series in Eqs. (34), here used 
to represent the Green functions, is assured by the general the- 
ory of the Green function (see, e.g., Smirnov 1964b) and is also 
proven in Appendix A. 



where the coefficient c„* are given by: 

= 2nF n f ° £ f(r,n)pr\\ - ^U^drdfi. 
Let us consider the equation: 

J (rf s) -Vf)dV = 

JVt \ nk ) 



(41) 



(42) 



(43) 
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obtained by means of Eq. (41). Considering Eqs. (30) and (38), 
Eq. (43) can be recast as: 



•JVr „ r n ,, 



(44) 



Moreover, if we introduce the volume average of the modulus 
of the perturbation |n| with respect to the weight function |V/|, 
i.e.: 



.JvX-'w 



then Eq. (44) gives a lower limit for it in the form: 



<|Ti C ' ,) |>V f 



■ Zk Cnkb n l 



(45) 



(46) 



The minimum dimensions of the volume Vf are set by the 
spatial resolution of the measurements of the angular velocity 
variations. They depend on the accuracy of the rotational split- 
ting coefficients, the inversion technique and the position within 
the convection zone (see, e.g., Schou et al. 1998; Howe et al. 
2005). The minimum order of the Jacobian polynomials N m 
needed to reproduce an angular velocity variation with a lati- 
tudinal resolution A8 is iV m ~ ||. Similarly, the minimum order 
of the radial eigenfunctions K m is set by the radial resolution Ar 
as K m ~ 2( ° R Ar rb) ■ Therefore, it is possible to truncate the series 
in Eqs. (44) and (46) to those upper limits for n and k because 
the coefficients c nk will decrease rapidly for n > N m and k > K m 
giving a negligible contribution to the sum. 

The statistical errors in the measurements of the angular 
velocity variations can be easily propagated through the linear 
equations (37), (39), (40) and (44) to find the errors on the esti- 
mates of V-Tj or the average of Ti . For instance, if we consider 
the standard deviations cr, of the data <f„ i.e., the rotational split- 
tings or the splitting coefficients from which the internal rotation 
is derived, the standard deviation crj of the integral in Eq. (44) 
is: 



(47) 



where 

e ink = 2nF n J ° f (l - ^p^c^^^drd/i, (48) 

and the functions c,(r,//) are the rotational inversion coefficients 
denned in Eq. (8) of Schou et al. (1998). 

Note that a constant relative error e = Aw/w in the measure- 
ments of A (e !) leads to the same relative error in Eq. (40) and in 
Eqs. (44) and (46), given the linear equations that relate the cor- 
responding quantities. As a matter of fact, there is also a system- 
atic error in our inversion method related to the poor knowledge 
of the turbulent viscosity %(r) that determines the form of the 
radial eigenfunctions £ nk . 



2.5. Kinetic energy variation and dissipation 

The variation of the kinetic energy of rotation associated with 
the torsional oscillations, averaged over the eleven-year cycle, 
can be computed after Lanza (2006b) and it is: 



(49) 



where 



{AT nk ) c = ^{[a^] 2 + [a ( :h 2 ) 



4F„ 



(50) 



The average dissipation rate of the kinetic energy of the torsional 
oscillations due to the turbulent viscosity is: 



n k 



(51) 



2.6. Torsional oscillations due to mean-field Lorentz force 

Most models of the torsional oscillations assume that they are 
due to the Lorentz force produced by the mean field as de- 
rived from dynamo models. Therefore, let us consider the case 
in which only the mean-field Maxwell stresses contribute to the 
perturbation, i.e.: 



1 

Ti = --(rsm6)B^B p . 
H 

If the mean radial B r and toroidal fields B$ are given by: 



B r = Z?o r cos(-crf), 



(52) 



Ba 



5cos(-<xf + n r ), 



(53) 



where LT r is the phase lag between the two field components, the 
components of T\ r in Eq. (31) are: 



' lr 



' lr 



1 r sin 6 

-- cos n r — - — BoiBo*, 

2 n 

1 . ,_, rsin# 

- sinlX — : — B Q[ B 04 ,. 

2 ii 



(54) 



An estimate of t^ c) and t\ s ' can be obtained from the method 
outlined in Sect. 2.4, considering a localization function f(r) that 
depends only on the radial co-ordinate. Specifically, Eq. (44) can 

(c) (s) 

be used to compute a volume average of t 1( . and t 1( . from which 
the average stress amplitude |B r Z?0| and phase lag IT can be de- 
termined. Analogous considerations are valid for the meridional 
component of the mean field B g and B^,. Adopting a localiza- 
tion function f(9) depending only on 6, it is possible to estimate 
IBeB^I and the phase lag Tig between Bg and B^. Such results are 
important to constrain mean-field dynamo models of the solar 
cycle, as discussed by, e. g., Schlichenmaier & Stix (1995) (see 
also Sect. 3.3). 



2.7. Application to solar-like stars 

Sequences of Doppler images can be used to measure the sur- 
face differential rotation of solar-like stars and its time variabil- 
ity, as done by, e.g., Donati et al. (2003) and Jeffers et al. (2007) 
in the cases of AB Dor and LQ Hya. Lanza (2006a) discussed 
the implications of the observed changes of the surface differen- 
tial rotation on the internal dynamics of their convection zones, 
assuming that the angular velocity is constant over cylindrical 
surfaces co-axial with the rotation axis. The present model al- 
lows us to relax the Taylor-Proudman constraint on the internal 
angular velocity, but some different assumptions must be intro- 
duced to obtain the internal torques in those active stars. Here we 
assume that the variation of the surface differential rotation is en- 
tirely due to the Maxwell stresses of the internal magnetic fields, 
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localized in the overshoot layer below the convection zone, as 
in the interface dynamo model by, e.g., Parker (1993). The in- 
terior model of the Sun can be applied also to AB Dor and LQ 
Hya because they have a similar relative depth of the convection 
zone. Therefore, the basic quantities can be scaled according to 
the stellar parameters, as explained in Lanza (2006a). Following 
Donati et al. (2003), we consider a surface differential rotation 
of the form: 



some average of the internal angular velocity of the star, say 
<jj m (f), can be measured as a function of the time: 



(O m (t) 



cu(r,/j., i)w(r, /u)drd/i, 



(60) 



Q(//, t) = Q e qW - dQ(t)//, 



(55) 



where w is an appropriate weight function that takes into account 
the averaging effects of the limited spatial resolution, and r\, is 
the radius at the base of the stellar convective envelope. Let us 
introduce an auxiliary weight function: 



where £2 eq is the equatorial angular velocity and dQ. the latitudi- w\ (r, jf) = 



nal shear, both regarded as time-dependent. Since P ( ^' l) = 1 and 

fi 2 = i + -j^j ^ can ^ e recast in the form of Eq. (16) leading 
to: 



3 k 



(56) 



where R is the radius of the star. The functions fi„k(t) can be 
computed by assuming that the timescale of variation of the dif- 
ferential rotation ?dr is significantly shorter than the timescales 
for angular momentum transport, as given by A~£. This as- 
sumption is justified in the case of rapidly rotating stars by the 
observed variation timescales of the order of a few years to- 
gether with the expected rotational quenching of the viscosity 
(see, e.g., Kichatinov et al. 1994). Therefore, Eq. (23) leads to 
a„k(t) - fDRj6«t(0- The angular velocity can be written in a form 
similar to Eq. (33) by introducing an appropriate Green function. 
For instance, considering the first of Eqs. (56), we find: 



i _ 3 r 

- -dQ = — 

5 8nJ v 



where: 

G s (r, r') = ^ £ok(r)£ok(r') 



(V-Ti)G s (fi,r'W, 



(57) 



(58) 



and the integration is extended over the stellar convection zone. 
Assuming that V • T\ is different from zero only in an overshoot 
layer of volume V a and applying considerations similar to those 
of Sects. 2.3 and 2.4, we find a lower limit for the variation of 
the averaged perturbation term: 



(61) 



which will not diverge toward the poles (fi = ±1) and the surface 
(p = 0) because the weight function w is localized into the stellar 
interior and goes to zero rapidly enough toward the rotation axis 
and the stellar surface. We can develop the function w\ as: 



AT™ K.,, 



(62) 



where the summation can be truncated at some low orders, say, 
N w and K w , because the function w\ is not sharply localized in 
the stellar interior. Considering Eq. (60), we find: 



(jj m (t) 



II 

n k 



1 



2nF„ 



WnkO!„k- 



(63) 



For the sake of simplicity, we assume now that the time scale 
of variation of the functions a„k(t) is long with respect to A~l so 
that Eq. (23) gives: (5„k - A n ka„k- Considering Eq. (30), we find: 



w m (f) 



I-lii 



V n 



W„k 

2nA nk 



(64) 



Equation (64) can be used to find a lower limit to the average of 
|ti| over the convection zone volume. If we indicate by M w the 
maximum of the function: 



w nk 



over the volume of the convection zone, we find: 
jylnldv | Wm(f )| 



<|Tl(f)l> = 



v 



(65) 



(66) 



<A(V-ti)> Vo > 



8tt 

3?DR 



AO, 



eq 



\A(dQ.) 



(59) 



where A£2 eq and A(c/Q) are the amplitudes of variation of the dif- 
ferential rotation parameters, and M is the maximum of G S (R, r') 
in the overshoot layer. 

This result made use of the limited information we can get 
from surface differential rotation. However, in the near future, 
the observations of the rotational splittings of stellar oscillations 
promise to give information on the internal rotation and its possi- 
ble time variations. Since only the modes of low degrees (I < 3) 
are detectable in disk-integrated measurements, the spatial reso- 
lution of the derived internal angular velocity profile is very low. 
Lochard et al. (2005) considered the case in which only the mean 
radial profile of £2 is measurable. 

In view of such an additional information accessible through 
asteroseismology, let us consider a more general case in which 



3. Application to the Sun 

3. 1 . Interior model, eigenvalues and eigenfunctions 

A model of the solar interior can be used to specify the func- 
tions p(r) and rj t (r) that appear in our equations. While the den- 
sity stratification can be determined with an accuracy better than 
0.5%, the turbulent dynamical viscosity is uncertain by at least 
one order of magnitude and it is estimated from the mixing- 
length theory according to the formula: 



1 



(67) 



where q'ml is the ratio of the mixing-length to the pressure scale 
height Hp and u Q is the convective velocity given by: 



/ a M LL Q y 
\40nr 2 pj 



(68) 
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0.70 0.75 0.80 0.85 0.90 0.95 

(r/Rs) 



Fig. 1. The ratio of the density to the density at the base of the 
convection zone (p/po- solid line) and the ratio of the turbulent 
dynamical viscosity to the turbulent viscosity at the base of the 
convection zone (r/,/rio - dotted) versus the fractionary radius 
r/R in our solar interior model. The density at the base of the 
convection zone is po = 0.1875 g cirT 3 and the turbulent viscos- 
ity is 770 = 2.56 x 10 12 g cirT 1 s -1 , respectively. A linear increase 
of % between the relative radii 0.673 and 0.713 is assumed to 
account for the effects of overshooting convection. 



where L Q is the luminosity of the Sun. In our computations we 
adopt a M L = 1 .5 and assume the solar model S for the interior 
quantities 1 (Christensen-Dalsgaard et al. 1996). In that model, 
the base of the convection zone is at r = 0.713/?©. We con- 
sider also the effect of an overshoot layer extending between 
r = 0.673/?© and the base of the convection zone within which 
the turbulent dynamical viscosity is assumed to increase linearly 
from zero up to the value at the base of the convection zone. The 
density and the turbulent viscosity are plotted in Fig. 1 where 
their values have been normalized at the values at the base of the 
convection zone, respectively. 

The basic equations of our model (i.e., 12, 13 and 14) can 
be made nondimensional by adopting as the unit of length the 
solar radius R G , as the unit of density po, i-e-, the density at the 
base of the solar convection zone, and as a unit of time to = 
poR^/rio, where 770 is the turbulent viscosity at the base of the 
convection zone. As a matter of fact, the value of the diffusion 
coefficients estimated from the mixing-length theory leads to a 
too short period for the solar cycle in mean-field dynamo models. 
Covas et al. (2004) adopted a turbulent magnetic diffusivity v t = 
3 x 10 11 cm 2 s _1 to get a sunspot cycle of ~ 1 1 yr. This implies 
rjo ~ pv t — 5.62 x 10 10 g cirT 1 s _1 in our model. 

The radial eigenfunctions and the Jacobian polyno- 
mials P^ 1 ' 1 ' have been computed from the respective Sturm- 
Liouville problem equations by means of the Fortran 77 subrou- 
tine sleign2.f 2 (Bailey et al. 2001). For the radial eigenfunctions 
(„k, the Sturm-Lioville problem has been solved with Neumann 
boundary conditions at both ends, set at 0.675/? o an d 0.99/?© to 
avoid divergence at the surface. For the Jacobian polynomials, 
limit point boundary conditions have been adopted at /i = ± 1 . 

Note that, from a rigorous point of view, it would be better to 
use the helioseismic estimate of at r = 0.99 where we fixed 

or 

our outer boundary for the computation of the radial eigenfunc- 
tions instead of the stress-free boundary condition that is valid 
only at the surface. However, the differences are confined to the 

1 See http://bigcat.ifa.au.dk/~jcd/solar_models/ 

2 http://www.math.niu.edu/SL2/ 



outermost layer of the solar convection zone, the moment of in- 
ertia of which is so small that there are no pratical consequences. 

The eigenvalues A n k and the eigenfunctions computed by 
sleign2.f have been compared with those computed by means of 
the code introduced by Lanza (2006b). The relative differences 
in the eigenvalues and in the eigenfunctions are lower than 1.5% 
for n < 14, k < 10. However, some problems of convergence 
of the numerical algorithm used by sleign2.f have been found 
for k > 20, particularly for n > 30, so we decided to limit its 
application up to k — 19. 

The Jacobian polynomials and the eigenvalues computed by 
sleign2.f are very good up to n = 30, as it has been found by 
comparison with their analytic expressions up to n = 10 and 
their asymptotic expressions for n > 12. We conclude that for 
n > 30 and k > 20 it is better to use the asymptotic formulae 
(22) and (18) instead of the numerically computed and P^' l) . 

The eigenvalue A n t gives the inverse of the characteristic 
timescale of angular momentum transfer of the mode corre- 
sponding to £„k under the action of the turbulent viscosity. The 
longest timescale corresponds to the lowest eigenvalue, i.e., 
/I20 = 0.078 in nondimensional units. It corresponds to a time 
scale of 0.86 yr with the turbulent viscosity given by the mixing- 
length theory and to 39.5 yr with 770 = 5.62 x 10 10 g cirT 1 s _1 . 

3.2. Localization functions for the source of the torsional 
oscillations 

The available data on the torsional oscillations are displayed 
with a typical radial resolution of 0.05 R Q and a latitudinal res- 
olution of 15° in Howe et al. (2005, 2006). The rotational inver- 
sion kernels of Schou et al. (1998) show a higher radial resolu- 
tion close to the surface, but, given the small amplitude of the 
torsional oscillations, the choice of a uniform resolution of 0.05 
R G seems to be better. 

We have found that the best results on the localization of the 
source term with the method outlined in Sect. 2.4 are obtained 
with localization functions that depend on r or /j only and that 
have a smooth derivative. As a typical function to probe the ra- 
dial localization, we adopt: 

10 for r < n, 

1 + sin [n(^)~ |] for r, < r < r 2 , (69) 
2 for r > Y2. 

This function has zero derivative, except in the interval Jn^t 
where its derivative varies smoothly reaching a maximum in the 
mid of the interval. In Fig. 2 we plot the modulus of the deriva- 
tive |^|, as given by the sum of the series of the radial eigenfunc- 
tions truncated at k = 19, for three different intervals centered at 
0.725, 0.85 and 0.965 R , respectively, all with an amplitude of 
r_ — rj = 0.05 R Q . The derivative is well approximated by the 
truncated series, with small sidelobes the amplitude of which in- 
creases toward the surface of the Sun because the eigenfunctions 
scale as (prj^^, according to the asymptotic expression (22). 

The performance of the localization method introduced in 
Sect. 2.4 has been tested with simulated data in the absence of 
noise. The results of two such tests are plotted in Fig. 3 where we 
consider the case of a purely radial perturbation T\ — r t r local- 
ized within an interval of amplitude 0.05 R e . Its time dependence 
is assumed to be purely cosinusoidal and the corresponding co- 
efficients (3„k are computed by means of Eq. (30) up to the orders 
11 = 38 and k = 19. From the /3 nJ t, the coefficients a„t are com- 
puted by means of Eq. (26) and the simulated angular velocity 
perturbation follows from Eq. (16). 
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0.90 0.95 



Fig. 2. The modulus of the radial localization kernel versus the 
relative radius for the function (69), as obtained by truncating 
the series of the radial eigenfunctions at K m = 19, for three dif- 
ferent intervals centered at r = 0.725 (upper panel), 0.85 (middle 
panel) and 0.965 R Q (lower panel), respectively, all with an am- 
plitude of r-i - r\ = 0.05 R Q . The modulus of the derivative has 
been normalized to its maximum value. 
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Fig. 3. Upper panel: Test case of the application of the inversion 
method introduced in Sect. 2.4. An input profile with T\ = 1.0? 
in nondimensional units between 0.80 and 0.85 R Q (plotted as 
the solid line) is used to simulate a noiseless profile of angular 
velocity perturbation. The reconstructed lower-limit profile ac- 
cording to Eq. (46) is plotted as a dotted line in the same panel. 
Lower panel: The same as in the upper panel, but with an input 
profile localized between 0.775 and 0.825 R Q . 



When the interval in which T\ is localized coincides with one 
of the inversion intervals [r\ , r2], the lower limit for \t\ | turns out 
to be ~ 80%-90% of the value assumed in the simulation (cf. 
Fig. 3 upper panel). The agreement increases up to 90%-95% 
if we consider a case in which n r has a constant sign and put 
the value of the derivative in the denominator of Eq. (46) 
instead of its modulus. This happens because the modulus of the 
derivative increases the effects of the sidelobes by increasing the 
value of the denominator in Eq. (46). When the interval in which 
the assumed T\ is localized does not coincide with an interval 
of the inversion grid, the inversion method still performs well 
distributing the contributions among neighbour intervals nearly 
in the correct proportion (cf. Fig. 3, lower panel). The case of 
simulations including a noise component is straightforward to 
treat thanks to the linear character of our inversion method. The 





0.0 



Fig. 4. The modulus of the latitudinal localization kernel versus 
p for the function (70), as obtained by truncating the series of 
the Jacobian polynomials at N m = 30, for three different inter- 
vals of p in the Northern hemisphere, i.e., [0, 0.26] (upper panel), 
[0.5,0.7] (middle panel) and [0.87,0.99] (lower panel), respec- 
tively. Note the symmetry of the modulus of the kernel with re- 
spect to the equator. The value of j ^ | is normalized at its maxi- 
mum at p = ± 1 . 

inverted value turns out to be the sum of the inverted noiseless 
value and of the contribution coming from the inversion of the 
noise. Their amplitude ratio is equal to the ratio of the amplitudes 
of the input noise to the input signal (for constant relative signal 
errors), as discussed in Sect. 2.4. 

The localization in latitude can be sampled by means of a 
localization function of the kind: 

(0 for -1 < p. < -p 2 , 

1+sin [*teH] f<* "AG <//<-//!, 
fip.) = {2 for-pi<p<pu (70) 

1+sin [!--(^r)] forp l <p<p 2 , 
for p2 < p < 1 . 

Such a function is symmetric with respect to the equator so that 
its derivative ^ is antisymmetric, as it is the source term T\ 
leading to a symmetric perturbation of the angular velocity. The 
localization function has a derivative always equal to zero ex- 
cept in two intervals ] - pi,-p\[ and ]p\,pi[, symmetric with 
respect to the solar equator. Its representation by means of the 
Jacobian polynomials with degree up to N = 30 is given in 

Fig. 4. Note that the maximum of j^j is reached in two sharp 

peaks at p = ±1. However, they are so narrow that their contri- 
bution to the integral in Eq. (46) is modest in comparison to those 
of the broader peaks in the intervals [-pi,-pi] and [p\,p2\. 
Several tests have been performed, as in the case of the radial 
localization, to assess the performance of the proposed method. 
The results are similar to those obtained for the radial case and 
are not discussed here. 

We conclude that our choice of N m = 30 and K m = 19 is 
perfectly adequate to invert the available data on the solar tor- 
sional oscillations to derive information on the location of the 
perturbation term within the convection zone. 

3.3. Results on the sources of the torsional oscillations 

The data plotted in Fig. 4 of Howe et al. (2006) can be used for 
an illustrative application of the inversion methods introduced 
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Fig. 5. The isocontours of the function A T (r, fj) as defined by 
Eq. (71) in the case of 770 = 2.56 x 10 12 g crrT 1 s _1 . The scale 
on the left indicates the ranges corresponding to the different 
colors and is in units of 10 5 g cirT 1 s~ 2 . The relative statistical 
uncertainty of A T is of ~ 30%, as it follows from the relative 
uncertainty of the data. 



in Sect. 2.4. They are given with a sampling of 15° between 
the equator and 60° of latitude, i.e., in the range in which the 
rotational inversion techniques perform better (cf. Schou et al. 
1998). The features distinguishable on the plots indicate an ac- 
tual radial resolution of ~ 0.05 R Q , in agreement with the sam- 
pling adopted in Fig. 3 of Howe et al. (2005). An average sta- 
tistical error of about 30% can be assumed for the amplitude, 
whereas the phase errors become very large below 0.80 R , es- 
pecially at low latitudes, because of the uncertainty in the re- 
construction of the signal in the deep layers. Note that the error 
intervals reported in Fig. 4 of Howe et al. (2006) indicate only 
how the particular method solution (here an OLA inversion of 
SoHO/MDI data) would vary with a different realization of the 
input data affected by a randon Gaussian noise. Unfortunately, 
they do not give the statistical ranges in which the true values 
of the amplitude and phase are likely to lie. This is not a major 
limitation in the context of the present study because we aim at 
illustrating the capabilities of the proposed approach rather than 
derive definitive conclusions. 

To perform our inversion, we interpolate linearly the values 
of the amplitude and phase over the grid used to compute the 
radial eigenfunctions and the Jacobian polynomials. We assume 
that the amplitude is zero at the poles and increases linearly to- 
ward 60° of latitude whereas the phase is constant poleward of 
60° of latitude. 

The divergence of the angular momentum flux perturbation 
Ti can be obtained from Eq. (40). We define its amplitude and 
phase as: 



a t = V[ v ' T i e) ] 2 + [ v ' T i 5) ] 2 ' 



<t> T = tan 



(V- 



r (l) , 



(71) 



(72) 



They are plotted in Figs. 



5 and 6, respectively, in the case of 
770 = 2.56 x 10 12 g crrT 1 s . A suitable smoothing has been ap- 



plied to have a resolution of ~ 0.05 R Q in the radial coordinate 
and « 15° in latitude. In order to show the effect of a smaller 
value of the turbulent viscosity, we plot in Figs. 7 and 8 the iso- 



0.67 ■ 




Fig. 6. The isocontours of <t> T (r, /i) as defined by Eq. (72) in the 
case of 770 = 2.56x 10 12 g cirT 1 s _1 . The phase ranges from — | to 
|. The scale on the left indicates the phase ranges corresponding 
to the different colors and is in radians. 
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Fig. 7. The same as Fig. 5 in the case of 770 = 5.62 x 10 g 
cm -1 s _1 . The scale on the left is in units of 10 3 g cirT 1 s~ 2 . The 
relative statistical uncertainty of A T is of about 30%. 
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contours of A T and O r for ?7 t = 5.62 x 10 g cm s 



The 



Fig. 8. The same as Fig. 6 in the case of 770 = 5.62 x 10 10 g cm 
s" 1 . 



amplitude of the perturbation term is higher close to the equator 
because the data in Fig. 4 of Howe et al. (2006) mainly sample 
the low-latitude branch of the torsional oscillations. The relative 
maxima of A T are reached close to the base of the convection 
zone and at a radius of ~ 0.85 R Q and their locations show only 
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a minor dependence on the value of 770. Conversely, the value 
of 770 significantly affects ® r , as it follows from Eq. (39). When 
r/o is large, the timescales for angular momentum exchange, as 
given by the inverse of the lowest eigenvalues, are significantly 
shorter than the eleven-year cycle and the perturbations are al- 
most in phase with the angular velocity variations. When 770 is 
sufficiently low, the timescales for angular momentum exchange 
become comparable or longer than the eleven-year cycle so the 
torsional oscillations lag behind the perturbations. Considering 
Fig. 6, we see that the phase is in agreement with that in Fig. 4 of 
Howe et al. (2006), except for r > 0.93 R e . The disagreement in 
those layers is due to the small values of V-Ti close to the surface 
that makes the corresponding phases uncertain (cf. Fig. 5). The 
phase lag is apparent by comparing Fig. 8 with Fig. 6, especially 
close to the equator for 0.72 < r/R Q < 0.92. It is less evident 
near the base of the convection zone and in the upper layers due 
to the smaller values of V • Ti in those regions. 

It is interesting to note that the dependence of the amplitude 
and phase of the torsional oscillations on the turbulent viscosity 
can lead to its estimate in the framework of mean-field mod- 
els (see, e. g., Riidiger et al. 1986; Riidiger 1989, for details). 
Specifically, Rempel (2007) finds that a mean turbulent kine- 
matic viscosity about one order of magnitude smaller than the 
mixing-length estimate is needed to reproduce the polar branch 
of the torsional oscillations. 

Lower limits for the modulus of the angular momentum flux 
vector |ti| can be derived from Eq. (46). From the lower limits 

on |t ( j C) | and |Tj' s) |, a lower limit on |ri | = ^[t^] 2 + [Tj S) ] 2 can be 
derived and it is plotted in Figs. 9 and 10 for the radial and lat- 
itudinal localizations described in Sect. 3.2, respectively. Given 
the uncertainty in the penetration depth of the torsional oscilla- 
tions, in addition to the results for a penetration down to the base 
of the convection zone, we plot also those for penetration depths 
of 0.80 and 0.90 R e , respectively. They are obtained simply by 
assuming that the oscillation amplitude has the value in Fig. 4 of 
Howe et al. (2006) above the penetration depth and drops to zero 
immediately below it. 

When the oscillations penetrate down to the base of the con- 
vection zone, the maximum of the perturbation term is reached 
between 0.80 and 0.85 R e . It is mainly localized into two latitude 
zones, i.e., within ±15° from the equator and between 30° and 
60°. Note that our data refer mainly to the low-latitude branch 
of the oscillations, so the possible source at latitude > 60° can- 
not be detected. When the turbulent dynamical viscosity is re- 
duced with respect to the mixing-length value, the amplitude of 
the perturbation term drops, but its radial and latitudinal local- 
izations are not greatly affected, except for a shift of the nearly 
equatorial band towards higher latitudes. 

Under the hypothesis that the Lorentz force is the only source 
of the torsional oscillations, the amplitude of the Maxwell stress 
B x B,p can be estimated from the lower limit of |ri r | and it turns 
out to be ~ 2.7 x 10 5 G 2 in the case of 770 = 2.56 x 10 12 g 
cm -1 s _1 . Considering that the poloidal field B r is of the order 
of 1 - 10 G, this leads to very high toroidal fields in the bulk of 
the convection zone, that would be highly unstable because of 
magnetic buoyancy. On the other hand, with 770 = 5.62 x 10 10 g 
cm" 1 s" 1 , the Maxwell stress is B^B^ 8.1 x 10 3 G 2 which leads 
to a toroidal field intensity of the order of 10 3 G. It may be sta- 
bly stored for timescales comparable to the solar cycle thanks to 
the effects of the downward turbulent pumping in the convection 
zone (cf., Brandenburg 2005, and references therein). 



Note that the maximum radial stress \B T B,p\ is reduced by a 
factor of ~ 30 when ?7 t is decreased from 2.56x 10 12 to 5.62X 10 10 
g cm" 1 s _1 , whereas the maximum of \BeB^\ is reduced only by 
a factor of ~ 4 (cf. Figs. 9 and 10). This mainly reflects the 
predominance of the radial gradient over the latitudinal gradient 
of the angular velocity perturbation in the deeper layers of the 
solar convection zone. 

The average phase lags n r and Tig between B r and Bg and 
the azimuthal field B^, as derived by the method in Sect. 2.6, are 
plotted in Figs. 11 and 12, respectively. It is interesting to note 
thatB r B > below - 0.85 R Q when 771 = 5.62xl0 10 g cm" 1 s"\ 
in agreement with the finding of most dynamo models in which 
the azimuthal field is produced by the stretching of the radial 
field in the low-latitude region, where > for r < 0.95 R e 
(cf., e.g., Riidiger & Hollerbach 2004; Schlichenmaier & Stix 
1995). Above ~ 0.85 R e , the phase relationship between the ra- 
dial and the azimuthal fields leads to B r B^ < with a phase lag 
of ~ n, in agreement with the early finding that the photospheric 
zone equatorward of the activity belt is rotating faster than that 
at higher latitudes (cf. Riidiger 1989). Note that ^ becomes 
negative for r > 0.95 R B at low latitudes, leading to a reversal 
of the phase relationship between the two field components in 
the outer layers of the solar convection zone. Our limited spatial 
resolution and the uncertainty of the measurements of the tor- 
sional oscillations inside the Sun might explain why we find the 
transition from a mostly positive to a negative B X B^, at ~ 0.85 
R e . 

The phase lag between B g and B^ depends remarkably on the 
colatitude for ?7 t = 2.56 x 10 12 g cm" 1 s"\ whereas it is almost 
constant for 77 t = 2.56 x 10 12 g cm" 1 s" 1 , leading in the latter 
case mostly to BgB$ < 0. This, together with ^ > 0, suggests 
that the toroidal field is mainly produced by the stretching of the 
radial field. 

On the other hand, if the angular momentum transport lead- 
ing to the torsional oscillations is produced only by a perturba- 
tion of the meridional flow, as in the thermal wind model, then a 
very small perturbation follows from the lower limit of |n|. For 
770 = 2.56 X 10 12 g cm -1 s" 1 , we find a minimum amplitude of 
the meridional flow component oscillating with the eleven-year 
cycle of « 3 cm s" 1 at a depth of 0.85 R Q and a latitude of 45°. 
Note, however, that if we estimate |n | from its divergence and 
the typical lengthscale of its variations in Fig. 5, we find a value 
about one order of magnitude larger, that is in agreement with 
the estimate by Rempel (2007). A similar argument applies also 
to the estimate of the Maxwell stresses made above. 

The average kinetic energy variation associated with the 
torsional oscillations, as computed from Eq. (49), is (7~) c = 
4.88xl0 28 erg, whereas the average dissipated power is 6.7xl0 26 
erg s" 1 , when 770 = 2.56 x 10 12 g cm" 1 s" 1 , and only 1.5 x 10 25 
erg s" 1 , when 770 = 5.6 x 10 10 g cm" 1 s" 1 . 

When the torsional oscillations are assumed to be confined to 
shallower and shallower layers, the location of the perturbation 
term is shifted closer and closer to the surface and it becomes 
more and more uniformly distributed in latitude. Its amplitude 
shows a remarkable decrease because of the smaller moment of 
inertia of the surface layers. 



4. Application to solar-like stars 

The results of Sect. 2.7 can be applied to the variation of the 
surface differential rotation observed in, e.g., LQ Hya between 
1996.99 and 2000.96 (see Table 2 of Donati et al. 2003). 
Assuming an internal structure analogous to the solar one and 
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„ 3.1 

'» 2.: 



77i=2.56 x 10 ,! g cm"' s" 



7)t-5.62 x 10 10 g cm-' s~ 



0.75 0.80 0.85 0.90 0.95 

(r/Ro) 

Fig. 9. Upper panel: The lower limit of the amplitude of the per- 
turbation term |Ti r | averaged over spherical shells of thickness 
0.05 R Q versus the relative radius; different linestyles and col- 
ors refer to the depth at which the torsional oscillations are as- 
sumed to vanish: black solid line - oscillations extending down 
to the base of the convection zone; green dotted line - oscilla- 
tions extending down to 0.80 R s ; red dashed line - oscillations 
extending down to 0.90 R e . Results are obtained assuming a tur- 
bulent dynamical viscosity at the base of the convection zone 
rjo = 2.56 x 10 12 g cirT 1 s _1 . Lower panel: As in the upper panel, 
but for 770 = 5.62 x 10 ll) g cirT 1 s _1 . The relative statistical un- 
certainties are in all the cases of about 30%, as it follows from 
the uncertainties of the data. 




Fig. 10. Upper panel: The lower limit of the amplitude of the 
perturbation term |tia| averaged over different latitude zones ver- 
sus their average value of p. Different linestyles are used to in- 
dicate the results for different penetration depths of the torsional 
oscillations, as in Fig. 9. Results are obtained assuming a tur- 
bulent dynamical viscosity at the base of the convection zone 



770 = 2.56 x 10 12 g cm 1 s Lower panel: As in the upper panel, 



but for 770 = 5.62 x 10 g cm s . The relative statistical un- 
certainties are in all the cases of about 30%, as it follows from 
the uncertainties of the data. 



that the differential rotation variations observed at the surface 
are representative of those at a depth of 0.99 R, we can es- 
timate a lower limit for |V ■ n| in the overshoot layer from 
Eq. (59). It is used to estimate the Maxwell stresses by assum- 
ing that: r -^p 6r\V ■ T\\, where r = 0.67 R is the radius 
of the overshoot layer and 6r = 0.04 R its thickness. In such 
a way, the minimum magnetic field strength turns out to be: 




T/ t =2.56 x 10 12 g cm" 1 s"' 



?7,=5.62 x 10'° g cm" 1 s" 



0.75 0.80 0.85 0.90 0.95 

(r/R s ) 

Fig. 11. Upper panel: The phase lag n r between B, and as 
derived by Eqs. (54) averaged over spherical shells of thickness 
0.05 R e versus the relative radius. Different linestyles are used 
to indicate the results for different penetration depths of the tor- 
sional oscillations, as in Fig. 9. Results are obtained assuming a 
turbulent dynamical viscosity at the base of the convection zone 
1 . Lower panel: As in the upper panel, 



770 = 2.56 xlO 12 gem 
but for 770 = 5.62 x 10 10 



gem 



T),=2.56 x 10 12 g err 



r,,=5.62 x 10'° g 1 



Fig. 12. Upper panel: The phase lag lie between Bg and B$ as 
derived by Eqs. (54) averaged over different latitude zones ver- 
sus their average value of p. Different linestyles are used to in- 
dicate the results for different penetration depths of the torsional 
oscillations, as in Fig. 9. Results are obtained assuming a tur- 
bulent dynamical viscosity at the base of the convection zone 
770 = 2.56 x 10 12 g cirT 1 s _1 . Lower panel: As in the upper panel, 
but for 770 



5.62 x K^gcm" 1 s _1 . 



#min = -ijBpBff, ~ 2500 G. However, if we assume a poloidal 
field strength of ~ 100 G, as indicated by the Zeeman Doppler 
imaging, we find an azimuthal field of ~ 6.2 x 10 4 G, which 
is in the range of the values estimated by Lanza (2006a) in the 
framework of the Taylor-Proudman hypothesis. Note that the re- 
sult is independent of 770 in the limit f DR <K A~} . Although this 
application is purely illustrative, it suggests that our method can 
be applied to derive estimates of the internal magnetic torques 
(as well as other sources of angular momentum transport) when 
asteroseismic results will become available, in combination with 
surface rotation measurements, to further constrain rotation vari- 
ations in solar-like stars. 
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5. Discussion and conclusions 



asymptotic limit, those series can be written as: 



We introduced a general solution of the angular momentum 
transport equation that takes into account the density stratifica- 
tion and the radial dependence of the turbulent viscosity r] t . Its 
main limitation is due to the uncertainty of the turbulent viscos- 
ity in stellar convection zones. It is interesting to note that the 
method of the separation of variables to solve Eq. (3) can be ap- 
plied also when rj t is the product of a function of the radius by 
one of the latitude. However, when % depends on the latitude, 
the angular eigenfunctions are no longer Jacobian polynomials. 

Our formalism can be applied to compute the response of a 
turbulent convection zone to prescribed time-dependent Lorentz 
force and meridional circulation. From the mathematical point 
of view, it is a generalization of that of Riidiger et al. (1986) and 
can be easily compared with it by considering that: P^'^Qx) = 
— 4 ^ - -^-P l n {0), where P n is the Legendre polynomial of 



. Moreover, our method can be used to 



«+2 d\i 

degree n and P\ - 
estimate the torques leading to the angular momentum redistri- 
bution within the solar (or stellar) convection zone, thus general- 
izing the approach suggested by Komm et al. (2003). The main 
limitation, in addition to the uncertainty of the turbulent viscos- 
ity, comes from the low resolution and the limited accuracy of 
the present data on solar torsional oscillations, particularly in 
the deeper layers of the convection zone. Actually, those lay- 
ers are the most important because torsional oscillations with an 
amplitude of ~ 0.5%- 1% of the solar angular velocity extend- 
ing down to the base of the convection zone lead to Maxwell 
stresses with an intensity of at least « 8 x 10 3 G 2 around 0.85 R Q 
or a perturbation of the order of several percents of the merid- 
ional flow speed at the same depth (cf., e.g., Rempel 2007). If 
the torsional oscillations are due solely to the Maxwell stresses 
associated with the mean field of the solar dynamo, we can es- 
timate also the phase relationships between B l7 B g and B^. Our 
preliminary results indicate that B^B^ > in the layers below 
~ 0.85 R Q and B^B^ < in the outer layers which, together with 
helioseismic measurements of the internal angular velocity, sug- 
gest that the toroidal field is mainly produced by the stretching 
of the poloidal field by the radial shear. 

Future helioseismic measurements may improve our knowl- 
edge of the torsional oscillations, essentially by extending the 
time series of the data or by means of space-borne instruments, 
like those foreseen for the Solar Dynamic Observatory (see, e.g., 
Howe et al. 2006). On the other hand, asteroseismic measure- 
ments may open the possibility of investigating similar phenom- 
ena in solar-like stars, particularly in those young, rapidly rotat- 
ing objects showing variations of the angular velocity one or two 
orders of magnitude larger than the Sun. 



nk A " k 
n 

where: 

h " s Z V). 



(A.l) 



(A.2) 



From the asymptotic formulae it follows that for a point 
in the domain [r^,R Q ] x [r\,,R Q ]x] - 1,1[, the quantity 

| yfF~„U(r)U(r')Pn 4 V)| is limited with an upper bound inde- 
pendent of n and k. Therefore, the series (A.2) that define the 
coefficients h„ are uniformly convergent in that domain if the 
series Y k y- converges. 

This can be proven by considering the inequalities (21) and 
the formula: 



1 U 



f-f ak 2 + y 2 \ y ) ' 

k=\ J 




(A.3) 



where z, a + and y are real numbers and the function g(x) is 
defined as: 



g(x) = n 



exp(7rx) + exp(-7rx) 



exp(7rx) - exp(-7rx) 



(A.4) 



Eq. (A.3) follows from the equality (1.217.1) of Gradshteyn & 
Ryzhik (1994). In this way, we find: 



2n(n + 3)g 

oo 

<y — 



n(n + 3)Ql 2 



k=\ 

M 



nk 



2n(n + 3)q 



n 2 P 



n(n + 3)ql 2 



n(n + 3)Ql 2 



n 2 P 



1 < 



(A.5) 



n 2 p 



n(n + 3)ql 2 
n 2 p 



This result indicates that h n 



: for n » 1 . Since the eigen- 



Vn(n+3) 

functions -sfF^P^'^ifi) form an orthonormal set, the conver- 
gence of the series in (A. 1) follows by the Riesz-Fisher Theorem 
given that the series Ym h„ ~ 2n converges. In such a way, 
the uniform convergence of the series in Eqs. (34) in the domain 
[r h ,R B ] x [r b ,R Q ]x] - 1, l[x] - 1, 1[ is proven. 
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Appendix A: Proof of convergence of the Green 
function series 

The convergence of the Green function series in Eqs. (34) can 
be studied by considering the asymptotic formulae for („k and 
the Jacobian polynomials, i.e., for n » 1 and k » 1. In the 
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